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Abstract 

The  purpose  of  this  contract  was  to  apply  various  mathematical  techniques  for  speeding  up  electro¬ 
magnetic  simulation  codes  to  present  day  commercial  codes.  The  techniques  to  be  investigated  in  the 
original  contract  were:  wavelet  methods,  the  non-uniform  fast  Fourier  transform,  and  multipole  meth¬ 
ods.  As  the  contract  evolved,  changes  were  made  in  two  ways:  first,  the  specific  techniques  investigated 
changed;  second,  the  contract  was  modified  to  incorporate  portions  of  another  MAFET  contract,  which 
was  part  of  the  larger  industrial  consortium.  That  contract  was  terminated  after  its  second  year  of 
funding  due  to  budget  constraints.  As  a  result,  this  contract  supported  the  remainder  of  that  contracts 
work.  Details  are  given  as  to  the  specific  work  performed,  goals  achieved,  and  future  work  to  be  carried 
out. 


1  Summary  of  Work  Performed 

The  work  carried  out  under  this  contract  is  broken  out  into  three  main  sections:  wavelet  methods,  improve¬ 
ment  of  matrix  condition  numbers,  and  fast  package  models.  Eachof  these  is  described  in  more  detail  in 
later  sections.  Brief  summaries  of  the  activities  are  now  given. 


1.1  Wavelet  Methods 

The  first  half  of  the  contract  (18  months)  concentrated  on  using  wavelet  methods  to  speed  up  the  matrix  solve 
of  moment  method  matrices.  Moment  methods  are  the  principal  technique  used  for  solving  electromagnetic 
planar  circuit  problems.  (The  other  major  technique  of  choice,  finite  elements,  is  more  popular  with  packaging 
problems.  It  is  examined  in  the  packaging  simulation  section  of  this  contract.)  The  investigation  concentrated 
on  getting  the  wavelet  methods  to  work  on  commercial  codes.  We  concentrated  on  HP-EEsof’s  Momentum 
product.  They  gave  us  a  copy  of  the  code  in  which  we  could  get  at  the  matrices  being  generated  in  the  solve 
portion  of  the  numerical  engine.  We  then  examined  converting  the  matrices  to  a  wavelet  basis,  and  solving 
the  resulting  matrix  in  the  wavelet  basis.  We  had  hoped  that  the  solve  would  be  extremely  fast,  owing  to  the 
matrix  being  made  more  sparse.  Unfortunately,  the  condition  number  of  the  matrices  were  poor.  A  poorly 
conditioned  matrix  cannot  be  solved  by  iterative  techniques  quickly,  as  the  solve  time  goes  as  the  condition 
number.  Indeed,  in  extreme  cases,  the  solver  will  never  converge,  when  the  finite  precision  arithmetic  of  a 
computer  is  included. 

We  therefore  had  to  develop  a  way  of  working  with  the  poorly  conditioned  matrices.  Fortunately,  the 
structure  of  the  matrices  after  they  have  been  converted  to  a  wavelet  basis  permits  a  fast  solve,  even  though 
they  are  badly  conditioned.  A  standard  solution  technique  for  matrices  is  the  LU  decomposition  method. 
Indeed,  this  is  the  standard  method  used  in  commercial  moment  method  codes.  If  the  matrix  has  N  entries, 
the  solve  time  goes  as  1 9(N3).  However,  the  method  will  terminate  in  a  finite  number  of  steps,  and  is 
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guaranteed  to  give  the  correct  answer.  (If  the  condition  number  is  poor  enough,  the  “correct”  answer  is 
meaningless,  due  to  roundoff  error.  A  direct  solve  helps  the  conditioning  problem,  but  doesn’t  get  rid  of 
it.)  We  developed  a  “fast”  LU  solve  for  a  matrix  which  has  been  converted  to  a  wavelet  basis.  It  works  the 
same  way  as  a  conventional  LU  solve,  except  we  used  the  structure  of  the  matrix  to  advantage  to  reduce 
the  solve  time  to  i9(N\nN).  The  method  is  therefore  as  fast  as  iterative  methods,  and  has  the  advantage 
of  working  better  on  ill-conditioned  matrices.  However,  it  must  be  emphasized  that  this  does  not  solve  the 
poor  conditioning  problem.  Accuracy  will  still  be  lost  due  to  round-off  errors.  Rather,  it  allows  us  to  work 
with  poorer  conditioned  matrices  than  we  otherwise  would.  This  work  is  summarized  in  the  paper  by  Gines, 
et  al.  [3] 

We  also  exploited  the  structure  of  the  wavelet  matrices  in  a  second  way.  The  matrices  created  by  moment 
methods  have  a  structure  which  depends  on  the  numbering  of  the  current  discretization  cells.  This  leads  to 
banding  of  the  matrices,  in  which  there  are  bands  and  subbands  of  larger  elements,  -and  near  zero  entries 
between  these  bands.  Conventional  matrix  solution  techniques  cannot  take  advantage  of  this  structure. 
Indeed,  the  banding  can  actually  be  deleterious  to  the  solve,  in  that  a  relatively  sparse  matrix  can  fill  up  in 
the  process  of  getting  the  solution.  We  therefore  developed  a  “two  dimensional”  wavelet  method.  The  idea 
is  that  the  wavelet  method  now  knows  about  the  banding.  The  matrix  remains  sparse  during  the  solve;  no 
fill  in  occurs.  The  method  is  summarized  in  Gines  et  al.  [2] 

1.2  Matrix  Conditioning  Methods 

The  matrices  that  are  produced  by  commercial  codes  are  typically  very  poorly  conditioned.  A  poor  condition 
number  makes  it  difficult  to  quickly  solve  the  matrix.  The  solve  time  of  iterative  methods  goes  as  the 
condition  number.  In  severe  cases,  it  is  not  even  possible  to  attain  a  solution  when  realistic  computer 
roundoff  error  is  considered.  Poorly  conditioned  matrices  will  also  hurt  wavelet  methods  for  matrix  solutions. 
A  wavelet  transformation  of  the  matrix  is  an  orthogonal  basis  transformation,  and  the  condition  number  will 
therefore  be  the  same.  This  in  turn  means  not  many  elements  of  the  wavelet  matrix  may  be  thrown  away, 
and  sparsification  of  the  matrix  will  not  be  that  significant.  Therefore,  little  speed  up  can  be  expected. 

We  examined  the  matrices  and  found  the  eigenvalue  structure  was  not  what  we  expected.  Normally,  a 
poorly  conditioned  matrix  will  have  eigenvalues  that  decrease  linearly  with  number,  until  the  final  eigenvalues 
axe  very  small.  The  condition  number  is  poor  as  it  is  the  ratio  of  the  largest  to  smallest  eigenvalues.  When 
we  looked  at  the  matrices  from  Momentum,  we  found  that  there  two  groups  of  eigenvalues.  The  first  group 
had  eigenvalues  that  were  all  about  the  same  order  of  magnitude;  the  second  group’s  eigenvalues  were  much 
smaller  than  the  first  group.  This  is  an  unusual  eigenvalue  structure  and  suggests  that  the  underlying 
equations  have  a  structure  that  can  be  exploited  to  improve  the  conditioning  of  the  matrix.  After  some 
searching,  we  found  that  other  researchers  had  similar  ideas.  In  particular,  we  were  influenced  by  the  work 
of  Vecchi  et  al  [7].  In  that  work,  Professor  Vecchi  explains  how  to  precondition  the  matrices  by  breaking 
them  up  into  two  parts,  and  treating  these  two  parts  differently.  The  original  matrix  is  converted  to  a  new 
basis  composed  of  “loops”  and  “stars”.  Loop  basis  functions  represent  a  current  which  has  no  divergence; 
star  basis  functions  have  no  curl.  By  breaking  the  current  up  in  this  manner,  it  is  possible  to  correlate  the 
loops  and  stars  with  the  two  groups  of  eigenvalues.  This  fact  can  be  exploited  to  solve  for  the  eigenvalues 
of  the  submatrices  (loop  -  loop  and  star  -  star).  The  eigenvectors  are  used  as  a  preconditioner,  whereby  the 
resulting  matrix  has  a  much  better  condition  number. 

Our  contribution  is  to  extend  some  of  Vecchi’s  ideas  to  work  more  practically  with  the  matrices  we 
typically  see  in  moment  method  codes.  In  particular,  Vecchi’s  method  relies  on  getting  a  static  solution  of 
the  submatrices.  We  get  the  solution  at  the  actual  frequency  of  interest.  We  also  believe  we  have  a  way  of 
getting  the  eigenvalue  structure  without  having  to  solve  the  complete  submatrix  problem.  This  work  is  still 
in  progress,  but  promises  to  significantly  increase  the  speed  of  the  method.  Details  are  given  in  the  technical 
section  of  this  report. 

1*3  Fast  Calculation  of  Parasitic  Coupling  in  Packages 

This  work  has  been  carried  over  from  the  other  terminated  MAFET  contract.  The  purpose  of  the  work 
is  to  be  able  to  include  parasitic  coupling  effects  in  circuit  simulators.  Parasitic  coupling  is  defined  to  be 
unintended  coupling  between  different  parts  of  a  circuit.  It  is  difficult  to  include  this  kind  of  coupling  in 
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a  circuit  simulator  because  it  is  usually  not  even  known  if  the  coupling  exists  or  what  it’s  level  is.  In 
earlier  work,  Dunn  et.  al  [4]  ,  developed  coupling  models  between  transmission  lines  in  a  circuit.  The 
model  consisted,  of  an  independent  voltage  and/or  current  source  in  the  coupled  line.  That  work  had  the 
limitation  of  needing  a  Green’s  function.  Therefore,  coupling  inside  a  package  cannot  be  calculated  except 
for  the  academic  exercise  of  an  empty,  rectangular  box.  This  work  extends  the  algorithm  to  include  arbitrary 
package  configurations.  This  is  accomplished  by  not  using  a  Green’s  function,  but  rather  relying  on  field 
data  from  a  finite  element  simulation  code.  In  a  sense,  the  finite  element  code  generates  a  numerical  Green’s 
function. 

We  are  still  working  on  results  for  the  method.  Details  are  given  in  the  technical  section  of  the  report. 

2  Publications,  Results,  and  Works  in  Progress 

This  section  list  the  publications  and  talks  resulting  from  this  work.  Work  in  progress  is  also  discussed.  No 
patents  resulted  in  this  contract. 

2.1  Wavelet  Methods 

The  principal  document  containing  the  wavelet  work  is  David  Gines  Ph.D.  thesis:  [2]  This  work  details  the 
direct  solve  method,  the  two-dimensional  wavelet  method,  and  a  host  of  examples. 

The  work  on  the  direct  solve  method  is  also  published  in  paper  form  as:  [3] 

The  work  on  two-dimensional  wavelets  is  being  reworked  in  a  more  computationally  desirable  form.  It  is 
expected  this  work  will  be  completed  this  academic  year,  at  which  point  a  paper  will  be  published. 

2.2  Matrix  Conditioning  Methods 

This  work  is  still  in  progress.  The  work  completed  to  date  has  been  presented  in  refereed  conference  publi¬ 
cations:  [6] 

We  will  be  presenting  the  latest  results  in  a  conference  this  spring  [5].  We  plan  on  submitting  a  full 
journal  paper  later  this  spring. 

2.3  Fast  Calculation  of  Parasitic  Coupling  in  Packages 

This  work  is  still  in  progress.  We  plan  on  submitting  a  paper  for  publication  later  this  spring. 

3  Technical  Description  of  Algorithms 

We  now  present  a  technical  description  of  the  algorithms  developed,  especially  for  the  work  where  publication 
is  in  progress. 


3.1  Wavelet  Methods 

Since  this  work  has  been  published  [3,  2]  ,  we  only  briefly  describe  it  here.  There  are  two  main  parts  to  this 
work:  fast,  direct  solution  of  waveletized  matrices;  and  a  two  dimensional  wavelet  algorithm. 

The  fast,  direct  solution  method  is  now  described.  A  full,  dense  matrix  is  converted  to  a  wavelet  basis 
in  the  standard  manner.  For  details  on  this  procedure,  see  for  example  [1],  In  particular,  we  used  the  non¬ 
standard  wavelet  transformation,  as  described  in  that  paper.  The  matrix  produced  has  the  standard  wavelet 
structure,  with  groupings  of  submatrices  following  the  diagonal.  The  matrices  we  obtained  from  Momentum 
are  typically  poorly  conditioned;  we  therefore  wanted  a  direct  method.  A  normal  LU  decomposition  on 
the  matrix  is  not  desirable  as  fill  in  will  occur  on  the  matrix,  which  is  mainly  zeros.  The  sparse,  direct 
method  developed  overcomes  this  method  by  recognizing  the  block  structure  of  the  matrix  and  preserving 
it  throughout  the  LU  decomposition  process.  Pivoting  is  implemented  within  the  blocks  to  keep  the  process 
conditionally  convergent;  thereby  rendering  it  useful  in  poorly  conditioned  matrices  examined  in  this  work. 
Details  are  found  in:  [2] 
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The  two  dimensional  wavelet  algorithm  again  works  on  the  waveletized  version  of  the  matrix.  Each  block 
of  the  matrix  has  subbanding  associated  with  the  numbering  of  the  original  cell  basis  of  the  problem.  This 
numbering  scheme  is  known,  and  therefore  the  subbanding  is  predictable.  The  two-dimensional  scheme  adds 
an  extra  set  of  indices  to  the  wavelets.  For  example,  we  have  four  indices  instead  of  the  normal  two  for  a 
one  spatial  dimensional  problem;  we  have  eight  indices  instead  of  four  for  a  two  spatial  dimensional  problem. 
This  process  is  accomplished  mathematically  by  working  with  matrices  of  matrices.  The  work  was  coded 
in  C.  The  code  is  complicated  because  of  the  complicated  indexing  required.  Benchmarks  showed  that  the 
method  will  be  faster  for  problem  sizes  above  5000  unknowns.  Below  that  size,  the  method  is  actually  slower, 
because  of  the  increased  overhead  required  in  the  indexing  scheme.  Details  are  given  in:  [2] 

3.2  Matrix  Conditioning  Methods 

The  goal  of  this  work  is  to  be  able  to  work  with  the  poorly  conditioned  matrices  resulting  from  the  moment 
method  solution.  First,  we  understood  the  reason  for  the  poor  conditioning.  We  used  the  work  of  Vecchi 
[8],  who  understood  the  essential  point.  The  integral  equation  used  is  an  electric  field  integral  equation 
(EFIE),  in  which  there  are  two  terms.  We  are  solving  for  the  unknown  current  on  the  signal  lines.  The 
first  term,  works  on  the  divergence  of  the  current,  or  equivalently  the  charge.  The  second  term  works  on 
the  whole  current.  The  integral  equation  can  have  solutions  which  are  nearly  in  the  null  space;  i.e.,  near 
zero  eigenvalues  are  possible.  In  particular,  these  eigenvalues  result,  because  there  are  currents  which  have 
zero  divergence.  A  zero  divergence  current  is  in  the  null  space  of  the  first  term;  the  second  term  gives  a 
contribution  proportional  to  the  frequency.  Therefore,  at  low  frequencies,  the  zero  divergence  current  is 
nearly  in  the  null  space  of  the  integral  equation  operators.  This  so  called  low  frequency  “catastrophe”  is 
evident  in  all  commercial  planar  solvers.  Unfortunately,  for  most  practical  problems,  we  are  close  to  this 
limit.  Even  if  a  solution  is  possible,  it  is  impossible  to  use  iterative  solvers,  and  sparsification  techniques,  for 
example  wavelet  methods,  will  do  little  good. 

Fortunately,  in  Vechhi’s  paper  a  method  is  suggested  to  fix  the  problem.  Specifically,  a  preconditioner 
is  developed  whereby  the  condition  number  of  the  transformed  matrix  is  lowered  substantially.  The  gist  of 
the  method  is  to  rewrite  the  current  into  two  different  sets  of  basis  functions:  divergence  free  and  curl  free. 
These  functions  are  usually  referred  to  as  loops  (divergence  free)  and  stars  (curl  free).  This  transformation 
on  the  matrix  is  useful  to  carry  out  because  of  the  two  terms  in  the  integral  equation.  The  first  term  only 
operates  on  the  loops;  the  second  term  operates  on  both  the  loops  and  the  stars.  Vecchi  showed  that  the 
resulting  matrix  has  a  very  clear  structure,  in  which  the  star-star  part  is  dominant,  and  the  loop-loop  part  is 
nearly  (but  not  totally)  independent.  This  suggests  a  preconditioner.  An  electrostatic  solve  of  the  problem 
is  made.  This  works  only  on  the  charge,  or  star-star  part  of  the  matrix.  This  problem  can  be  solved  using 
fast,  electrostatic  methods.  The  resulting  eigenfunctions  can  be  used  as  a  new  basis  for  the  star  part  of  the 
current.  Only  a  few  are  needed.  The  loop-loop  part  of  the  problem  is  preconditioned  in  the  same  way.  A 
magnetostatic  solve  is  carried  out,  with  the  resulting  eigenfunctions  being  used  as  the  preconditioning  basis. 
Vecchi  has  demonstrated  that  the  method  works  well  for  patch  antenna  problems. 

Vecchi’s  method  has  some  drawbacks  for  our  applications.  First  of  all,  he  requires  a  static  solve  of  the 
problem,  which  is  awkward  when  we  are  trying  to  use  solutions  from  commercial  codes.  Second,  the  static 
solves  can  take  t?(7V3)  amount  of  time.  The  method  is  therefore  expensive  to  run.  It  does  bring  down  the 
condition  number,  but  only  results  in  time  savings  if  several  frequency  points  are  used,  and  the  same  basis 
can  be  reused.  We  therefore  set  out  to  solve  these  two  drawbacks.  The  first  problem  is  overcome  by  taking 
the  star-star  and  loop-loop  blocks  of  the  matrix  and  performing  a  singular  value  decomposition  on  them 
(SVD).  This  method  allows  us  to  get  preconditioning  bases  at  the  frequencies  of  interest  without  having  to 
use  a  static  solution.  We  then  reuse  the  basis  over  the  frequency  band  of  interest.  Benchmarks  show  that  the 
method  is  useful  in  reducing  the  condition  number  over  the  frequency  band  of  interest,  with  one  caveat.  It 
is  difficult  to  go  through  a  resonance  in  the  frequency  response  of  the  circuit.  Details  are  given  in:  [5].  This 
solution  fixes  with  static  solve  issue,  but  we  still  have  speed  problems.  SVD  is  a  very  expensive  algorithm: 
tf(5N3).  Fortunately,  we  only  have  to  run  it  at  one  frequency  in  our  band.  We  have  a  new  idea  how  to  better 
solve  this  problem.  We  are  investigating  using  QR  algorithm,  instead  of  SVD.  The  advantage  of  this  is  that 
the  QR  algorithm  gives  us  the  most  “important”  eigenvalues  and  eigenvectors  first.  We  don’t  need  to  get 
all  the  eigenvalues.  Recall  that  the  final  basis  has  only  a  few  dominant  terms.  We  can  get  these  terms  first, 
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and  stop  the  algorithm  quickly.  This  is  the  same  idea  as  is  used  in  some  image  compression  techniques.  The 
question  is  whether  QR  is  stable  for  the  matrices  we  are  using.  We  think  it  will  be,  but  we  still  have  to  test 
the  procedure.  We  hope  to  have  definitive  results  by  later  this  spring. 

3.3  Fast  Calculation  of  Parasitic  Coupling  in  Packages 

The  algorithm  is  based  on  previous  work  by  Dunn  [4].  A  method  is  described  in  that  paper  where  the 
coupling  between  two  arbitrary  oriented  lines  on  a  substrate  is  calculated.  The  algorithm  works  by  exciting 
one  of  the  lines,  and  observing  the  fields  it  produces  at  the  other  line.  These  fields  produce  voltages  and 
currents  in  the  second  line.  These  same  voltages  and  currents  can  be  reproduced  by  modeling  the  coupling  as. 
independent  voltage  and  current  sources  in  the  second  line.  The  model  is  approximate  in  that  no  interaction 
back  to  the  first  line  is  accounted  for,  and  energy  is  not  conserved  in  the  sense  that  the  independent  voltage 
and  current  sources  create  energy.  Nonetheless,  it  can  be  a  useful  approximation  if  the  coupling  is  relatively 
weak,  is  as  typically  the  case  in  parasitic  coupling  problems.  The  independent  voltage  and  current  sources 
can  be  placed  in  a  circuit  simulation  package,  and  give  the  designer  an  idea  if  coupling  is  a  problem. 

The  main  drawback  of  this  work  is  that  it  requires  a  Green’s  function  in  order  to  calculate  the  fields  at 
the  second  line.  Therefore,  generalized  packages  are  not  allowed.  We  have  generalized  the  method  by  using 
a  finite  element  program  (HP-HFSS)  to  calculate  the  fields  from  the  first  line.  The  fields  can  be  extracted 
anywhere  in  the  space  of  the  simulation,  in  particular  where  the  second  line  is.  The  independent  voltage  and 
current  sources  are  then  calculated.  An  integration  is  carried  out  over  the  second  line,  where  the  magnetic 
or  electric  field  from  the  first  line  is  in  the  integrand.  The  method  is  limited  in  the  accuracy  obtained  for 
the  fields.  Typically,  large  tetrahedra  are  used  in  a  region  where  the  second  line  exists.  We  are  still  testing 
the  method  on  a  variety  of  canonical  structures.  Preliminary  data  show  that  the  method  works  well  enough 
that  order  of  magnitude  estimates  of  coupling  problems  can  be  made.  Notice,  that  the  advantage  of  this 
method  over  having  the  finite  element  program  simply  calculate  the  coupling  between  two  lines,  is  that  now 
the  second  line  can  be  anywhere  in  the  package.  Only  one  simulation  need  be  carried  out  per  package.  This 
is  in  contrast  to  requiring  a  new  simulation  every  time  the  second  line  is  moved.  We  expect  to  be  publishing 
the  results  later  this  spring. 


References 

[1]  G.  Beylkin,  R.  R.  Coifman,  and  V.  Rohklin.  Fast  wavelet  algorithms  and  numerical  algorithms  I.  Com- 
munications  in  Pure  and  Applied  Math ,  44:141-183,  1991. 

[2]  David  L.  Gines.  Fast  Electromagnetic  Simulations  Using  Wavelets .  PhD  thesis,  University  of  Colorado, 
Boulder,  May  1997. 

[3]  David  L.  Gines,  John  M.  Dunn,  and  Gregory  Beylkin.  LU  factorization  of  non-standard  forms  and  direct 
multiresolution  solvers.  Journal  of  Applied  Computational  and  Harmonic  Analysis ,  5:156—201,  1998. 

[4]  Lincoln  C.  Howard,  Kent  Larson,  and  John  M.  Dunn.  An  efficient  numerical  algorithm  for  the  calculation 
of  coupling  between  lines  in  MICs.  IEEE  Transactions  on  Microwave  Theory  and  Techniques ,  41(8):1287- 
1293,  August  1993. 

[5]  Hugh  MacMillan  and  John  M.  Dunn.  A  practical  method  for  increasing  the  speed  and  stability  of  the 
matrix  solve  in  moment  method  simulations  within  a  frequency  band.  Applied  Computational  Electro¬ 
magnetics  Conference  (ACES),  Monterey,  CA ,  March  1999. 

[6]  Hugh  MacMillan  and  John  M.  Dunn.  Eigenvalue  studies  of  matrices  resulting  from  EFIE  simulations  for 

planar  structures.  Applied  Computational  Electromagnetics  Conference  (ACES),  Monterey,  CA ,  March 
98.  ■ 

[7]  Giuseppe  Vecchi  et  al.  A  numerical  regularizatoin  of  the  efie  for  three-dimensional  planar  structures.  Int 
J  Microwave  Millimeter-Wave  CAE ,  7:1-22,  February  1997. 


5 


[8]  Giuseppe  Vecchi  et  al.  Analysis  of  cavity  modes  for  the  efie.  IEEE  Trans,  on  Antennas  and  Propagation, 
AP-44: 234-239,  April  1998. 


6 


